Variational ground states of the two-dimensional 
Hubbard model 



D Baeriswyl, D Eichenberger and M Menteshashvili 

Department of Physics, University of Fribourg, CH-1700 Fribourg, Switzerland. 
E-mail: dionys . baeriswylOunif r . ch 

Abstract. 

Recent refinements of analytical and numerical methods have improved our 
understanding of the ground-state phase diagram of the two-dimensional (2D) Hubbard 
model. Here we focus on variational approaches, but comparisons with both Quantum 
Cluster and Gaussian Monte Carlo methods are also made. Our own ansatz 
leads to an antiferromagnetic ground state at half filling with a slightly reduced 
staggered order parameter (as compared to simple mean-field theory) . Away from half 
filling, we find d-wave superconductivity, but confined to densities where the Fermi 
surface passes through the antiferromagnetic zone boundary (if hopping between both 
nearest- neighbour and next-nearest-ncighbour sites is considered). Our results agree 
surprisingly well with recent numerical studies using the Quantum Cluster method. 
An interesting trend is found by comparing gap parameters A (antiferromagnetic or 
superconducting) obtained with different variational wave functions. A varies by an 
order of magnitude and thus cannot be taken as a characteristic energy scale. In 
contrast, the order parameter is much less sensitive to the degree of sophistication of 
the variational schemes, at least at and near half filling. 
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1. Introduction 

The conventional mechanism of Cooper pairing rests on phonon exchange, but the case of 
superfluid 3 He shows that no additional degrees of freedom are necessary for establishing 
pairing, even in a system of fermions with predominantly repulsive interactions. Several 
years before the discovery of superfluidity in 3 He Kohn and Luttinger used perturbation 
theory to show that even a purely repulsive bare interaction can lead to an effective 
attraction in a channel with large enough angular momentum [I]. An analogous situation 
may prevail in superconducting cuprates, as suggested by Anderson as early as 1987 
[2]. The Hubbard Hamiltonian with a large on-site repulsion U is indeed the natural 
model for describing the antiferromagnetic Mott insulator of undoped cuprates. Clearly 
there are phonons in these materials and the electron-phonon coupling can be rather 
strong, as suggested by optical spectroscopy [3], but a thorough treatment of both 
strong correlations and electron-phonon coupling is required to understand to what 
extent phonons affect the electronic properties in the cuprates. A recent debate [H [5] 
on the origin of a kink observed in photoemission experiments demonstrates that such an 
analysis is still lacking. Correspondingly, the role of phonons in the superconductivity 
of cuprates remains also unclear. 

In this paper we restrict ourselves to the two-dimensional (2D) Hubbard model, 
described by the Hamiltonian 

H = H + UD, (1) 

where 

#o = -^VU> ( 2 ) 

ijcr 

describes hopping over the sites of a square lattice and 

D = s ^n^n ii (3) 

i 

measures the number of doubly occupied sites. The operator c ia (ctj annihilates 
(creates) an electron at site i with spin a and n ia = c\ a Ci a . Only hopping between 
nearest (tij = t) and next-nearest neighbours (t^- = t') will be considered. 

For small U the approach of Kohn and Luttinger looks promising, but naive 
perturbation theory diverges and one has to sum entire diagram classes (see 
[6] for a pedagogical treatment of this problem). Moreover, there are several 
competing instabilities close to half filling, in particular, rf-wave superconductivity 
and antiferromagnetism, or rather spin-density waves. The problem can be solved in 
an elegant way using the functional renormalization group [3 El El [10], which treats 
the competing density-wave and superconducting instabilities on the same footing. 
The specific techniques used by the different groups to calculate effective vertices and 
susceptibilities differ in details, but the overall results are more or less consistent, with 
an antiferromagnetic instability at half filling and c?-wave superconductivity away - but 
not too far away - from half filling. 
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In the large-t/ limit, where double occupancy is suppressed, the Hubbard model 
can be replaced by the Heisenberg model at half filling [IT] and by the t — J model close 
to half filling [12], defined by the Hamiltonian 

Ht-j = -J2 tijPo4 a c ja P + J v ( S i ■ S j ~ \ n i n J J > ( 4 ) 

ija ij 

where rij = + counts the number of electrons at site i, the projector 

p =n( i - n ^) (5) 

i 

eliminates configurations with doubly occupied sites, Si are spin | operators and 

Jij = mu. 

While the ground state of the 2D Heisenberg model is fairly well understood [13] - 
it is widely accepted that it exhibits long-range antiferromagnetic order with a reduced 
moment due to quantum fluctuations - the ground state of the t — J model remains 
an unsolved problem, despite an extensive use of sophisticated methods during the last 
two decades [HJ [15] . A simple ansatz 

|*) = A|*o>, (6) 

is frequently used, where \^q) is the ground state of a suitable single-particle 
Hamiltonian, representing for instance the filled Fermi sea, a spin-density wave or a 
superconductor. In a large region of doping the most favourable mean-field state \^/ ) has 
been found to be a superconducting ground state with rf-wave symmetry [T6 l fT71 fT8l fT9] . 
In this case the ansatz ([6]) describes a resonating valence bond (RVB) [20] state, jokingly 
referred to as the 'plain vanilla version of RVB' |21j . 

Sometimes the t — J model is used to provide a simple argument for pairing. If the 
exchange coupling J is larger than the hopping amplitude t then it is favourable for two 
holes to remain nearest neighbours. Unfortunately, this is an unphysical limit for the 
Hubbard model for which the t — J model is only a valid approximation for J <^t. 

Experiments on layered cuprates indicate that U is neither small enough for 
a perturbative treatment nor large enough for the mapping to the t — J model. 
Thus, neutron scattering experiments on undoped cuprates [22j can be very well 
interpreted in terms of the 2D Hubbard model with U between 8 1 and 10 1 (depending 
on the approximation used for the magnon spectrum [23j |21]). in agreement both 
with an analysis of Raman data [25] and with a comparison between theoretical and 
experimental optical absorption spectra [26J. As to the 'kinetic' term H , angular- 
resolved photoemission experiments [27l [28] give evidence for a hole-like Fermi surface, 
which can be described by taking both nearest-neighbour and next-nearest-neighbour 
hopping terms into account [29], [30]; the experiments are consistent with t! m —0.3t. 

A lot of effort has been spent to treat the Hubbard model in the crossover regime 
of moderately large U. Quantum Monte Carlo methods have been developed in the 
1980s, but they suffer from severe statistical uncertainties at low temperatures and for 
large system sizes, both for many-fermion and quantum spin Hamiltonians [31]. This 
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'minus sign problem' turns out to be intrinsically hard [32]. The recently introduced 
Gaussian Quantum Monte Carlo method [33] leads to positive weights, but the method 
is not immune to systematic errors [HI]. Quantum Cluster methods, which replace the 
single site of the Dynamical Mean-Field Theory [35] by a cluster of several sites, offer a 
promising alternative route for the intermediate-?/ regime [36], [37]. However, the spatial 
extent of the cluster that is treated essentially exactly has so far been very small (2x2). 
The large changes found in the case of the Mott transition when proceeding from a 
single site to a (2 x 2) cluster [38] indicate that larger cluster sizes are needed to obtain 
accurate results. 

Alongside Quantum Monte Carlo and Quantum Cluster methods, variational 
schemes have played a major role in the search of possible ground states of many- 
body systems. In the context of the Hubbard model, many variational wave functions 
go back to the ansatz of Gutzwiller [39] 




where the variational parameter g reduces double occupancy, weakly for small U and 
completely for U — > oo. In spite of the simplicity of the ansatz, a substantial part of 
the correlation energy is obtained for small values of U. At half filling, the variational 
ground-state energy tends to the exact limiting value for U — > oo, where double 
occupancy is suppressed (g — > oo) and the ansatz (JTj) coincides with the RVB state (jSJ). 
Nonetheless many properties such as spin-spin correlation functions are described rather 
poorly, especially in the large-f/ regime. Several improvements have been proposed 
in terms of additional operators in front of the Gutzwiller ansatz, either to enhance 
magnetic or charge correlations [4"0l |4"T] or to increase the tendency of empty and doubly 
occupied sites to be next to each other [121 03]. Very recently more elaborate states 
have been proposed, one including backflow [H], the other introducing a large number 
of variational parameters, both into |\&o) and into the (Gutzwiller- Jastrow) correlation 
factor [32]. 

The main purpose of this paper is to describe our own variational ansatz in some 
more detail than previously [4"6l Wf\ and to compare our results with those of other 
approaches. In section [2] our ansatz is defined and shown to recover Anderson's RVB 
state in the large-t/ limit and for small hole doping. Section [3] analyses various wave 
functions for the Hubbard model on a four-site plaquette. For the square lattice, a 
variational Monte Carlo method is used to study our ansatz. The procedure is explained 
in section H] and applied to an antiferromagnetic state at half filling. Our results for 
superconductivity away from half filling are presented in section [5] and compared to 
those of other variational wave functions and also to a very recent computation using 
a Quantum Cluster method. We do find evidence for a superconducting ground state, 
with an interesting hint at a magnetic mechanism. The concluding section [6] gives a 
brief summary of our main results. 




(7) 
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2. Variational ground states 

The primary effect of the Gutzwiller variational state (EJ) is the increased suppression 
of double occupancy as a function of U. Alternatively, we may motivate the ansatz as 
follows. Consider a general Hamiltonian H = Hq + H- m t, where Hq is a single-particle 
operator and H mt the interaction term. The state 

|*a> = e-^ltft) (8) 

tends to the exact ground state of H for A — > oo (or to one of the ground states 
in the case of degeneracy) unless the trial state \^ t ) is orthogonal to |\&). If \^ t ) does 
not differ too much from a small parameter A may be sufficient to obtain a good 
approximation. In this case we can make the replacement 

e" A " « e" A ^ e"^ . (9) 

Choosing |\P t ) = l^o); the ground state of i^o; an d treating A as a variational parameter, 
we recover the Gutzwiller ansatz (0) in the case of the Hubbard model. At the same 
time this derivation indicates that the Gutzwiller ansatz is best suited for small values 
of U. 

The factorization ([9]) is not unique, we may also choose 

e -XH ^ e -XH e -XH int _ ( 1Q j 

With |^ t ) — |^o) we arrive at an ansatz where both exponentials have to be kept. 
Allowing them to vary independently, we obtain 

|*GB> = e-* t e-^|* >. (11) 
The operator e~ 9D partially suppresses double occupancy for g > 0, while e~ hH °^ 
promotes both hole motion and kinetic exchange. The limit h — > leads back to 
the Gutzwiller ansatz. The variational state ( TTTi) has been introduced by Otsuka [IE] 
and studied both in one and in two dimensions. He found a substantial improvement 
with respect to the Gutzwiller ansatz. Moreover, |^gb) has a large overlap with the 
exact ground state for all values of U in one dimension, whereas on a square lattice this 
is only true for relatively small values of U (smaller than the bandwidth). 

For g — > oo and h <C 1 the ansatz (TIT!) leads to the t — J model. To see this, 
we consider the limit U — > oo, where e~ 9D \^o) is replaced by the r.h.s. of ([6]). The 
expectation value of the energy is then given by 

E _ (fr |P e-"W f (ff + Ufye-^PolVo) 
^ \P e-^/tp \^ ) 

At half filling the parameter h is equal to —t/U in the large-?/ limit [19] and vanishes 
for U — > oo. This remains valid very close to half filling, as we readily demonstrate. 
Expanding the expression ( fT2l) in powers of h, we obtain for the numerator 

(^olAe-^tfo + UD^-^P^v) = $! Q \P Q H Q P Q \V Q ) - 2^(V \P H 2 P \V ) + 

+ U^(^o\PoH PiH P \^ ) + - , (13) 
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where Pi is a projector onto the subspace with one doubly occupied site. For small 
doping the operator Hq in the second term can be replaced by H PiH , and the 
denominator in ( Tl"2l) by 1. The minimization of the energy with respect to h then 
gives h = —t/U and 

Po|*o>- (14) 
Close to half filling, the second term in this expression can be rewritten as 

PoffoAffoA « 2 4 " St ■ SA . (15) 

hi 

Thus we find indeed that in the limit U — > oo the variational treatment of the Hubbard 
model using the ansatz pip is equivalent to the variational treatment of the t — J model 
using the ansatz ([6]). However, we have to keep in mind that the use of (I lip for g — > oo 
is of doubtful validity Therefore variational results obtained with the fully projected 
Gutzwiller wave function, as in the 'plain vanilla' RVB theory, have to be handled with 
care. 

A complementary approach [19] starts from the ground state |^oo) for U — > oo, 
ideally the ground state of the t — J model close to half filling. Applying the same 
arguments as above, we obtain two new variational states, 

|*b> =e-^°/ t |tt 00 ), 

|*BG) = e- fl6 e- ft ^/ t |* 0O ) l (16) 

which should be well suited for describing the large- C/ regime. Unfortunately, it is in 
general very difficult to deal with the variational states \^b) an d |^bg) ; because even 
an approximate calculation of |^oo) represents already a highly nontrivial problem [50J. 

The trial states (JTj), ( TTTT) and ( jl~6l) yield an appealing picture for the Mott metal- 
insulator transition as a function of U at half filling [51], [52J, [531 [51] , where both |\&g) and 
|^gb) represent metallic states, while \^b) and |^bg) are insulating. For the soluble ID 
Hubbard model with long-range hopping, where tij ~ \i — the variational energies 
for |^ G ) and |\&b) are dual to each other at half filling; similarly |^gb) and |^bg) form 
a dual pair [5Tj . On the metallic side (U smaller than the bandwidth) |^gb) yields an 
order of magnitude improvement for the ground state energy as compared to |^g); and 
the same effect is observed on the insulating side when \^b) is replaced by |^bg)- 

3. Lessons from the Hubbard square 

The study of small clusters can provide useful benchmarks for variational wave functions. 
Clearly the cluster size should not be too small. In fact, for two particles on two sites 
both |^g) and \^b) yield the exact ground state, and there is no need to use more 
elaborate wave functions. Therefore we proceed to the Hubbard square, an essentially 
ID system when hopping is restricted to pairs of adjacent sites (figure[Ha)). We consider 
first the case of four particles and denote the corresponding ground state by |\I/( 4 )). 
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Figure 1. (a) Illustration of the Hubbard square; (b)-(d) representations of the three 
different types of configurations present in the ground state of the half-filled Hubbard 

square. The symbols stand for: O empty site, • • singlet bond, fi doubly 

occupied site. 



A theorem by Lieb [55] directly implies that the ground state of the ID Hubbard 
model for iV particles on iV sites is non-degenerate for any positive value of U and 
has total spin 5 = 0. For U — the ground state of our Hubbard square is fourfold 
degenerate (for S z = 0) with energy E = —At, and one can use perturbation theory to 
determine the U — > limit of |\p( 4 )). Apart from S = 0, |^^) can be characterized 
by its transformation properties with respect to the symmetry group T>^ it has total 
quasimomentum ir (i.e. it is odd under the rotation of the square by the angle 7r/2 
about the axis perpendicular to the plane of the square and passing through its centre), 
it is even with respect to reflections in a plane parallel to one of the sides and odd with 



respect to reflections in a plane passing through a diagonal. Denoting d\ :- 



with 

|$2> = 

l$o>: 



1 

V2 
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and 



bji, the ground state can be written as 



a[$2> + + c|$ ) 
Va 2 + b 2 T^ 



(17) 
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where |0) is the vacuum (zero-electron) state and a,b,c are real, [/-dependent 
coefficients, whose analytical expressions are too cumbersome to be presented here. The 
states |$d), illustrated in figure [TJb)-(d), are eigenstates of D with D = 0,1,2 doubly 
occupied sites (doublons). With increasing U the weights of and |$ 2 ) decrease and 
vanish for U — > 00. Therefore |$o) is equal to |^oo), the normalized ground state of 
the Heisenberg model on the square with nearest-neighbour antiferromagnetic exchange. 
We also note that in the state |$ 2 ) the two doublons sit on adjacent sites of the square. 
Similarly, in the state |$i) the doublon and the holon (empty site) sit on adjacent sites 



Variational ground states of the 2D Hubbard model 



8 



2 particles 



0.001 




0.5 1 

U/4t 



0.5 

4t/U 



0.5 

4t/U 



Figure 2. Relative deviation AE/E of different variational ground-state energies from 
the exact value E for the Hubbard square: (a) half filling, (b) quarter filling. The lines 
correspond to: |^ G ), |* B >, |*gb)- 



of the square. The hopping operator Hq mixes the states \$d), but does not bring in 
other states. 

As already noted by other authors [56], the ground state l^/^) has an affinity for d- 
wave pairing. In fact, the quantity (^ (4) |C d C]|^ (4) ), where C\ := | (&! 2 — &23 + & L — & 4i) 
creates a <i-wave pair, is larger than the corresponding quantity for s-wave pairing, 
(* (4) |C s Cl|* (4) ) with C\ := |(&i2 + & 23 + 44 + & 4i)- Tne difference, especially prominent 
for small U, decreases with increasing U and tends to zero for U oo. The ground 
state |\l/( 4 )) also exhibits local antiferromagnetic order, which is noticeably enhanced 
with increasing U. 

Knowing the exact ground state, one can examine the quality of various variational 
wave functions. Figure[2](a) shows that the Gutzwiller ansatz \^g) already yields a fairly 
good approximation for all values of U, with the biggest deviation for intermediate U 
(~ 4t), where the variational energy deviates from the exact value by about 0.5%, 
whereas the overlap (\1/g|^^ 4 ' ) ) is still about 99.95%. The state I^b), equation ffl6|) . is 
even better, with an improvement of one order of magnitude (or more) for U > It 
as compared to |^g)- Generally speaking, since the wave function ( ITTjl has two 
independent coefficients, one needs in principle a variational wave function with at least 
two parameters to reproduce the exact ground state. Not all wave functions are equally 
successful, though. Both the refined Gutzwiller ansatz |\&gb) and its counterpart |^bg) 
do give the exact ground state |\I/( 4 )) of the half-filled Hubbard square, but Kaplan's 
generalization [12] fails in doing so, because what should be suppressed by his additional 
operator (configurations with doublons having no holons on nearest-neighbour sites) is 
absent for the states |$d), D = 0, 1, 2. More involved variational states [H] can easily 
reproduce the exact ground state. 

The ground state \^^) for two electrons on the four sites of a square is rather 
different. Its energy reaches the limiting value = —2y/2t (and not 0) for U — > oo. 
Its quasimomentum is and the two particles form an s-wave pair. The quantity 
(\|/( 2 )|(7 s (7t|v|/(2)^ expressing the affinity for s-wave pairing, dominates now over the 
corresponding <i-wave expression, (^W\C d C\\^W), for all values of U except for U = 0, 
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where the two quantities are equal. The difference increases with U and reaches a 
maximum for U — > oo. 

I 1 !^ 2 )) can also be represented as a superposition of three components, but now only 
one of them contains real-space configurations with doubly occupied sites. As a result, 
the Gutzwiller ansatz \^g) gives a fairly poor description of the state, it even predicts a 
wrong asymptotic value for the energy, E G — « 0.16 1 for U — > oo (figure EJ^b)). The 
refined Gutzwiller ansatz |^gb) is n °t perfect either, it reproduces the exact ground 
state only below a limiting value of U « 3.46 1. Even if one adds a third operator e~ 9 ' D 
in front of I^gb); the situation remains essentially the same. What needs to be done in 
order to achieve full coincidence with the exact ground state for all U is to use PqHqPq 
instead of H in combination with the Gutzwiller factor, thereby ensuring that no more 
doublons are created by H after being suppressed by e~ 9D . 

On the other hand, |^b) gives a more or less satisfactory approximation of the exact 
ground state, about as good as the Gutzwiller ansatz in the half-filled case. Adding 
to I^b) any of the above-mentioned factors (Gutzwiller, doublon-holon correlations, 
Jastrow) gives full coincidence with the exact two-electron ground state. 

The analysis of the Hubbard square confirms that proceeding from |^g) to |^gb) 
(or from |^b) to |^bg)) improves substantially the variational ground state. Moreover, 
for U > t it is preferable to use |^b) (I^bg)) as a variational ansatz, rather than |^g) 
(I^gb)), i-G. it is advantageous to start from |\l/oo)> the ground state of the t — J model 
in the limit J — > (or of the Heisenberg Hamiltonian at half filling). 



4. Antiferromagnetism 

We return now to the square lattice. The lessons from the Hubbard square suggest the 
use of |^b) or I^bg); but, unfortunately, the reference state |^oo) is very difficult to 
handle, except for small clusters (using exact diagonalization). This is the main reason 
why |^b) an d I^bg) have not been used so far for the square lattice. In the remainder 
of this paper we will therefore mostly be concerned with wave functions linked to \^o), 
the ground state of an appropriate mean-field Hamiltonian. 

We first consider an antiferromagnetic ground state, restricting ourselves to the 
half-filled band case (number of particles N equal to the number of sites L). To enforce 
a commensurate spin-density wave (with a preference for up spins on one sublattice and 
for down spins on the other), we introduce the mean- field Hamiltonian 

H m{ = H - A ]T(-l)> lT - nil ) , (18) 

i 

where % is even on one sublattice and odd on the other. The staggered order parameter 
m is defined as 

m:= ^E(- 1 ) l ^T-^l), (19) 

i 

where the expectation value is calculated with respect to the chosen trial state. The 
simplest case is the ground state of H mi , a commensurate spin-density wave (SDW) with 
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wave vector Q = (jr, ir). 

The expectation value of the Hubbard Hamiltonian can be written as 

(H) = (H mf ) + 2LAm + U ^(n*^) ■ ( 20 ) 

i 

For the trial state |SDW), the ground state of H m {, the average double occupancy is 
simply given by 

(SDW|n iT n a |SDW) = \- m\ (21) 
while the Hellman-Feynman theorem yields 

— (SDW|iJ mf |SDW) = -2Lm. (22) 

Minimizing (H) with respect to A and using ( 1211) and ( |22i) . we obtain 

Htt? 

(A-C/m) — = 0, (23) 

and therefore 

A = Um. (24) 

The order parameter m is easily calculated for the mean-field state. The 
transformation 

Cfccr = 7= / G ^ Qct j (25) 

V ^""^ 

where the wave vectors k belong to the first Brillouin zone and Ri, i = l...,L, are lattice 
vectors (with lattice constant set to 1), diagonalizes the hopping term, 

H = ^2e k c{ a c ka , (26) 
with a spectrum 

£fc = — 2i(cos k x + cos ky) — At' cos k x cos k y . (27) 

In the rest of this section, we choose t' — and thus obtain a spectrum with electron-hole 
symmetry, 

£fc = —£k±Q ■ (28) 
In reciprocal space, the mean-field Hamiltonian reads 

Hwt = Yl' { £ k(cl a c ka - c ] k+Qa c k+Q ^) - ( jA(c t fc(7 c fc+ Q iCT + c{ +Q(T c kcT )} , (29) 

fccr 

where the notation ^' fc means that the sum includes only half of the wave vectors of 
the Brillouin zone, those satisfying e k < 0. H mi is easily diagonalized by a Bogoliubov 
transformation. The spectrum (in the folded zone) is split into ±E k , where 



E k = J el + A2 . (30) 
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The order parameter is found to be 

k K 

The relations (1211) and f[3"Tj) determine both the gap parameter A and the order 
parameter m for the mean-field ground state |SDW). 
We now proceed to the correlated trial state 

|*gb) = e~ hko,t e~ g£) I SDW) . (32) 

To decouple the terms n^n^ in the operator e~ 9D , a discrete Hubbard-Stratonovich 
transformation [57] is applied, 



3 ~9 D = exp <J -g ^ n H n n 



JJ cosh [2a (n iT - ra a ) ] exp | - 1 (ra iT + n a ) J 
2~ L exp <^ ^(2ao-ri - |)n; CT I , 

ri,...,r £ I, iff J 



(33) 



where r, are Ising variables assuming the values ±1, and a = arctgy / th((yf/4). As a result, 
both exponential operators in ( 132]) are quadratic in fermionic creation and annihilation 
operators, and therefore the fermionic degrees of freedom can be integrated out [58J. 
We obtain the following expressions 

(* gb |£|*gb}= Yl E(t u ...t l ;t[,...,t' l ), 



■n.,—,TL 

r[,...S L 



(*gb|*gb> = Yl N(r 1 ,...r L -y i ,...y L ), (34) 



T 1 ,...,T L 



where the quantities E(ti, ...Tl', t[, r^) and N(ti, ...t l ;t[, ...,t' l ) are products of 
determinants of L/2 x L/2 matrices. (Initially, one has to deal with L x L matrices, 
but these can be block-diagonalized into two L/2 x L/2 matrices because up and down 
spins are never mixed.) 

The expectation value of the energy 



(^GB 


#|^GB> 


(^GB 


^gb) 



E(g,h,A)= \r 7 (35) 



is then calculated by Monte Carlo sampling and minimized with respect to the three 
variational parameters g, h, A. The computations have been limited to an 8 x 8 lattice, 
with a fixed value of U = St. 

It is instructive to follow the evolution of the results as the variational ground state 
is refined, from the simple spin-density wave state \^>q) = |SDW) via the Gutzwiller 
ansatz \^ G ) = e"^|SDW) to |tf GB ) = e _ft ^°/*e- flA |SDW). The gap parameter A is 
found to decrease dramatically [46]. It amounts to about 3.6 1 for the plain SDW state, 
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1.3 1 for the Gutzwiller ansatz and 0.32 1 for | \&gb) • This indicates that A has no clear 
physical meaning and cannot be simply related to an excitation gap for the correlated 
states |^g) and |\I/gb)- We expect A to decrease steadily as the variational ansatz 
is further refined and to tend to zero as the trial state approaches the exact ground 
state. In fact, in this limit the gap parameter merely plays the role of an (infinitesimal) 
symmetry-breaking field. 

In spite of the strong variation of A, the order parameter does not change much |46j . 
from m pa 0.45 for |SDW) tom« 0.43 for |^g) an d m ~ 0.39 for |\1/gb)- Thus the result 
for our most sophisticated variational ansatz is higher than the extrapolation m pa 0.22 
from Monte Carlo simulations on relatively small lattices [59] , and also higher than the 
value m ~ 0.31 extracted from large-scale Monte Carlo simulations [60] for the 2D spin- 
| Heisenberg model (and expected to represent an upper limit for the antiferromagnetic 
order parameter of the Hubbard model). 

Quite generally, the order parameter obtained within mean-field theory is reduced 
when quantum fluctuations are taken into account. Our results indicate that quantum 
fluctuations are progressively included when proceeding from the simple spin-density 
wave state |SDW) over \^g) to |\Pgb)- Nonetheless, long-range (spin wave) fluctuations 
are suppressed by the gap A, which is small but still finite for |\Pgb)- It is interesting to 
note that a recent calculation using a Quantum Cluster method [61] gave m pa 0.4, which 
agrees with our value. The small size of the cluster (2 x 2) used in these calculations 
suggests that long-range fluctuations are not properly taken into account either. 

5. Superconductivity 

We discuss now the region away from half filling, i.e. densities n = N/L ^ 1. We have 
tried to examine the fate of antiferromagnetism for n ^ 1, but so far our variational 
Monte Carlo approach did not converge fast enough to allow the extraction of reliable 
results. In the following we limit ourselves to superconducting ground states with d-wave 
symmetry. The mean-field state \^q) is constructed as a conventional Bardeen-Cooper- 
Schrieffer (BCS) state, with 



The variational calculations follow the same steps as in the SDW case, but now one has 
to deal with 2L x 2L matrices because, in contrast to the SDW case, up and down spins 
are now mixed as are particles and holes. 
Our trial state 



has three variational parameters, h, g and A , as well as the parameter /i, which fixes 
the average number of electrons, but is not identical to the true chemical potential. 
Instead of this 'grand-canonical' set-up one could also work with a BCS state projected 
onto a fixed number of particles, where /i becomes a fourth variational parameter. 
Unfortunately, the minus sign problem turns out to be severe in the 'canonical' case 



A(fc) = A (cos k x — cos k y ) . 



(36) 



^ GB ) = e-^V^lBCS) 



(37) 
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[52] , presumably because |BCS)tv is a correlated state, whereas the conventional BCS 
wave function can be written as a single Slater determinant. The results discussed below 
have all been obtained using the grand-canonical version. 

The calculations were mostly done for an 8 x 8 lattice, but in order to study the 
size dependence we have also made a few runs for 6x6 and 10 x 10 lattices. The results 
for an 8 x 8 and a 10 x 10 lattice do not differ considerably if the gap parameter A 
is large enough [16], as is the case for n = 0.90 and 0.94. For n = 0.84 the finite size 
effects are still rather strong. 

As in the case of antiferromagnetism, the gap parameter A varies strongly as the 
trial wave function is refined. Thus, the fully Gutzwiller-projected BCS state, used as 
an ansatz for the t — J model [18] , gives large values in an extended region of densities. 
Within this approximation, the largest gap parameter, Ao ~ t, is found at half filling, 
followed by a nearly linear decrease as a function of doping concentration x = 1 — n, 
until A vanishes at x c ~ 0.35. Similar behaviour is found for \^g) applied to the 
Hubbard model [46], but with a smaller critical density and a reduced value at x — 0, 
A ~ 0.2 1. The gap is reduced still further when proceeding from \^g) to |\I/gb) 
[46] , with a maximum of 0.13 1 for x ~ 0.1 and a critical hole density x c ~ 0.18. It 
is worthwhile to mention that a Gutzwiller ansatz with coexisting antiferromagnetism 
and superconductivity also produces a maximum in the gap parameter as a function of 
density [63], A ~ 0.1 1 at x « 0.1. 

For the repulsive Hubbard model, the superconducting order parameter $ is 
commonly chosen as the expectation value of a pair of creation operators on neighbouring 
sites, 

$ = (4 T cJ +ri ) . (38) 

For d-wave symmetry, $ has a different sign for a horizontal bond than for a vertical 
bond. Alternatively, one can extract the order parameter from the correlation function 

S TT <(Ri - Rj) = {c i+Ti Ci^c ] j+Tl[ ) (39) 

for \Ri — Rj\ — > oo. For the exact ground state, the two procedures, the evaluation of 
( 1381) for an infinitesimal symmetry breaking on the one hand, and the determination 
of the asymptotic behaviour of the correlation function (|39|) on the other, are 
expected to yield the same order parameter $ in the thermodynamic limit (N, L — > 
oo for constant density n). 

We first discuss results for if = 0. Figure [3] shows the <i-wave order parameter 
as a function of doping, x = 1 — n, obtained following four different routes. Three 
of them are variational, full Gutzwiller projection for the t — J model (black squares) 
|18j . the Gutzwiller ansatz, including a finite antiferromagnetic order parameter close 
to half filling (red triangles) [63J, and our ansatz fl37|) (green circles) [46]. The null 
result (blue diamonds) is from recent (Gaussian) Monte Carlo simulations of Aimi and 
Imada [64], who extracted tiny upper bounds for the order parameter from the long- 
distance behaviour of the correlation function (|39|) . Unfortunately, the four data sets 
were obtained for four different values of U . Nonetheless, the three variational results 




Figure 3. Comparison of superconducting order parameters for the simple Hubbard 
model and four different approaches: full Gutzwillcr projection for the t — J model 
(U = 12 1, ■ — ■), partial Gutzwillcr projection (U = 10 1, ▲ — A), our ansatz (U = 
8t, • • ), and Gaussian Monte Carlo (U = 6t, ♦). 



are remarkably similar for weak doping (x < 0.1). For x > 0.18 there appears to be a 
discrepancy between the results of [63] and our results. Note that the Gutzwiller data 
[63] have been obtained on the basis of the correlation function ( |39i) . which is expected 
to yield a finite order parameter for a finite system, even if there is no long-range order in 
the thermodynamic limit. Therefore the apparent discrepancy may be an artefact of the 
procedure, as pointed out in [63]. This is confirmed by our own results for the Gutzwiller 
trial state \^g)i where we find [58] that the extrapolation of the order parameter for 
L — > oo gives $ — > for x > 0.2, in agreement with our results for |^gb)- Clearly 
our ansatz |^gb) should be superior to the Gutzwiller wave function, either partially or 
fully projected. It is worthwhile to mention that refined variational states for the t — J 
model [50] did not give markedly different results from those presented in figure [3] for the 
fully Gutzwiller-projected wave function. Therefore there seems to be a clear difference 
between the Hubbard and t — J model predictions for superconductivity above, say, 
a doping concentration of 18%, where superconductivity appears to be absent for the 
Hubbard model, but to persist up to about 40% for the t — J model. This is not a true 
discrepancy because the mapping from the large- i7 Hubbard model to the t — J model 
is only justified close to half filling. Figure [3] also shows that there is no disagreement 
between our variational result [16] and the Gaussian Monte Carlo data [63], which so 
far are only available for x > 0.18. 

We turn now to the more realistic case of hopping between both nearest and next- 
nearest neighbours. Figure H] compares our variational results (red circles) [37] with 
recent data obtained with the Quantum Cluster method (open triangles) [61] . The same 
parameter sets (U = 8t, t' — — 0.3 1) were used in both approaches. The agreement 
is remarkably good, in view of the fact that the two methods are very different. Only 
the surprising peak found in the Quantum Cluster approach for an electron doping of 




Figure 4. Superconducting order parameter as a function of doping for the Hubbard 
model including next-nearest-ncighbour hopping: • — • our variational data, A — 
A results obtained with the Quantum Cluster method. All data points were obtained 
using the parameter values U = 8 1, t' — —0.3 1. 



about 5% is absent in our variational data. This peak is suppressed if antiferromagnetic 
long-range order is taken into account |61j . 

As noted previously [37], the interval —0.2 < x < 0.25, to which superconductivity 
is confined according to our variational study, corresponds to densities where the Fermi 
surface passes through the boundary of the folded Brillouin zone, at so-called hot 
spots. This means that parts of the Fermi surface are connected by wave vectors close 
to Q — (vr, 7r), for which the magnetic structure factor has a pronounced peak [37] , 
reminiscent of long-range antiferromagnetic order at half filling. This observation lends 
support to a magnetic pairing mechanism. In fact, an effective attraction of the form 
— (U 2 /t) Xo{k + k') is deduced using second-order perturbation theory for the effective 
particle-particle vertex [10], where Xo(<?) is the spin susceptibility for non-interacting 
electrons and the incoming and outgoing electrons have wave vectors ±k and ±fe', 
respectively. A similar expression for the effective attraction, —\{U /t) x(q, w), has been 
fitted to the particle-particle vertex calculated with the Quantum Cluster method [65J, 
with a density- and temperature-dependent coupling constant U and the numerically 
calculated spin susceptibility x(q, 



6. Summary and concluding remarks 

In this paper, we have analysed various variational wave functions in the context of 
the 2D Hubbard model. Both antiferromagnetic ground states (at half filling) and 
superconducting ground states with (i-wave symmetry (away from half filling) have been 
studied. We have considered states of the general form 

|*) = P|*o>, (40) 
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where |\&o) is the ground state of an appropriate mean- field Hamiltonian, which includes 
a variational parameter A, the gap parameter. The operator P is unity in the simple 
mean-field theory, it consists of a single operator, P = e~ 9D , in the Gutzwiller ansatz or 
of several operators in more refined wave functions. The optimized gap parameter was 
found to depend sensitively on the choice of the wave function. It decreases by about an 
order of magnitude when proceeding from the simple mean-field theory to our favourite 
ansatz with P = e- flH o/ t e -9D _ Therefore A, taken as a variational parameter in wave 
functions of the form ( l40i) . cannot be identified with a physical energy gap (except in the 
simple BCS case where A is the minimum energy for adding an electron to the system), 
nor can it be associated with a characteristic energy scale such as k^T* (where T* is the 
so-called pseudogap temperature [ZD US]). In contrast, the order parameter was found 
to depend little on the sophistication of the variational ansatz. 

To our knowledge, there is no rigorous proof for the existence of long-range 
antiferromagnetic order in the 2D Hubbard model at half filling (nor in the case of 
the spin~2 Heisenberg model), although this question is not strongly debated. The 
case of (d-wave) superconductivity is much more controversial. Both the relatively 
small disparities between the order parameters obtained with different variational states 
close to half filling and the generally good agreement between our results and those of 
completely different approaches, Quantum Cluster [61] and Gaussian Monte Carlo [64|, 
give rather strong support for the existence of (i-wave superconductivity in the Hubbard 
model on a square lattice. Nevertheless, as we have seen in the simple example of the 
Hubbard square, our wave function is probably not optimal for the large value of U 
used in our calculations (and believed to be appropriate for the cuprates). Progress 
with complementary wave functions (linked to the U — > oo limit) or with other methods 
would be highly desirable. 
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